#include "SQ_ZZ_pEXFactoring.h"

#include <fstream>

#include "tools.h"

#include <NTL/mat_ZZ_pE.h>
#include <NTL/ZZ_pXFactoring.h>


/** 
 * Given a set of integers with specific multiplicity, decomposes an integer in terms of the elements of the set
 *
 * e.g. given s=12 and the set of pairs {{6,1},{3,4},{1,1}}, where {6,1} means there is one 6, returns [1,2,0] as in 12 = 6*(1) + 3*(2) + 1*(0) 
 *
 * @return 1 if found a solution, 0 otherwise
 *
 **/


int int_decomposition(int* i,int* e_i,int s,int num_of_ints,int*& solution)
{
  // base conditiona
  if( s == 0 )
    {
      solution[num_of_ints-1]=0; // mark the current place with a 0
     // if(num_of_ints >=1 )
//	int_decomposition(i,e_i,0,num_of_ints-1,solution); // call the function again to automatically 0 the rest
      return 1;
    }
  if(num_of_ints == 1)
    {
      for(int j=0;j<=e_i[num_of_ints-1];j++)
	{
	  if(s-j*i[num_of_ints-1] == 0) // found a solution
	    {
	      solution[num_of_ints-1] = j; // record the proper multiplicity
	      return 1;
	    }
	  else if( s-j*i[num_of_ints-1] < 0)
	    break;
	}
      return 0;
    }

  int j=1;
  for(;j<=e_i[num_of_ints-1];j++)
    {
      if(s-j*i[num_of_ints-1] >= 0 )
	{
	  int status =  int_decomposition(i,e_i,s-j*i[num_of_ints-1],num_of_ints-1,solution);
	  if( status == 1 ) // found a solution
	    {
	      solution[num_of_ints-1] = j;
	      return 1;
	    }
	}
      else // this element does not contribute to finding a solution
	{
	  solution[num_of_ints-1] = 0;
	  break;
	}
    }
  return  int_decomposition(i,e_i,s,num_of_ints-1,solution);
}



/**
 *	evaluates a polynomial in x at another polynomial, i.e. g(b) mod f
 **/
ZZ_pEX	SQ_eval(const ZZ_pEX& f,const ZZ_pX& g,const ZZ_pEX& b)
{
	ZZ_pEX res,Q;
	for(int i=0;i<=deg(g);i++)
	{
		if( !IsZero(coeff(g,i)) )
		{
			ZZ_pEX temp;
			SQ_divr(SQ_power(b,i),f,Q,temp); 		
//cout << "temp: " << temp << endl;
			res+=coeff(g,i)*temp;
		}
	}	
//cout << "res: " << res << endl;
	return res;
}

long MYgenerate_GF4(vec_ZZ_pE& element)
{	
	ZZ_pE alpha;

	char str1[256],str2[256];	// inefficient way of generating the primitive element
	ifstream In_tempfile;
	In_tempfile.open("primitive_element.ntl");
        In_tempfile.getline(str1,256);
	In_tempfile.close();
//cout << "given: " << str1 << endl;
//fflush(stdout);

	while(1)
	{
	         alpha = random_ZZ_pE();
		 if( !IsZero(alpha) && !IsOne(alpha) ) 
		{
			ofstream tempfile;
			tempfile.open("primitive.element");
			tempfile << alpha << endl;
			tempfile.close();

			ifstream tempfile1;
			tempfile1.open("primitive.element");
			tempfile1.getline(str2,256);
			tempfile1.close();			
//cout << "str2: " << str2 << endl;
			if(!strcmp(str1,str2))
				break;
		}
	}
//cout << "alpha: " << alpha <<  endl;

	long num_of_elements = to_long(ZZ_pE::cardinality());
	element.SetLength(num_of_elements);
	clear(element[0]); // element[0] = 0
	set(element[1]); // element[1] = 1
	element[2] = alpha;
	for(int i=3; i<num_of_elements;i++)
		element[i]=element[i-1]*element[2];

	return num_of_elements;
}


/**
 * given a row-reduced echelon form, solves the homogenous system of equations Ax=0 with GF(2) as fixed field
 */
/*
vec_ZZ_p OLD_solve_equ_GF2(mat_ZZ_pE A,int* leading_ones_pos)
{
	// solution vector
	vec_ZZ_p b; // b=[0 ... 0] is the soluction vector as in Ax=b
	int m = A.NumRows();
	int n = A.NumCols();
	b.SetLength(n);
	ZZ_p p_ONE;	// p_ONE = 1
	set(p_ONE);
	int mask[n];	// mask to keep track of changes
	memset(mask,-1,sizeof(int)*n);
	
	for(int i=m-1;i>=0;i--)
	{
		if( !IsZero(A[i]) )
		{
			mask[leading_ones_pos[i]]=1;
			for(int j=n-1;j>leading_ones_pos[i];j--)
			{
				if( !IsZero(A[i][j]) )
				{
					if( mask[j]==-1 ) 
					{
						b[j] = p_ONE;
						mask[j]=1;
					}
					b[leading_ones_pos[i]] += p_ONE;		
				}
			}		
		}
cout << "mask: " << endl;
for(int i=0;i<n;i++)
	cout << mask[i] << " ";
cout << endl;
cout << "b: " << endl << b << endl;
	}
	for(int i=0;i<n;i++)	// set independent variables to 1
		if(mask[i]==-1)
			b[i] = p_ONE;

	return b;
}
*/

/**
 * given a row-reduced echelon form, solves the homogenous system of equations Ax=0 with GF(2) as fixed field
 */
vec_ZZ_p solve_equ_GF2(mat_ZZ_pE A,int* leading_ones_pos)
{
	// solution vector
	vec_ZZ_p b; 
	int m = A.NumRows();
	int n = A.NumCols();
	b.SetLength(n);

	ZZ_p p_ONE;	// p_ONE = 1
	set(p_ONE);

	int mask[n];	// mask to keep track of changes
	memset(mask,-1,sizeof(int)*n);

	vec_ZZ_pE element;
	long order_of_GFq = MYgenerate_GF4(element);


	int available_alpha_pos[n];
	int available_1_pos[n];
	int available_alpha2_pos[n];
	int count_1,count_alpha,count_alpha2;
	
	for(int i=m-1;i>=0;i--)
	{
		count_1 = count_alpha = count_alpha2 = 0;

		if( !IsZero(A[i]) )
		{
			mask[leading_ones_pos[i]]=1;	// leading ones have to be assigned, they are not independent
			ZZ_pE rowb; // b=0 is the soluction vector as in Ax=b for each row

			if(i==0)
				rowb = p_ONE;
				

			for(int j=leading_ones_pos[i]+1;j<n;j++)
			{
				if( !IsZero(A[i][j]) )
				{
					if( mask[j]==-1 ) 
					{
						if(A[i][j] == element[1])	// == 1
						{
							available_1_pos[count_1] = j;
							count_1++;
						}
						else if(A[i][j] == element[2]) // == alpha
						{
							available_alpha_pos[count_alpha] = j;
							count_1++;
						}
						else // == alpha^2
						{
							available_alpha2_pos[count_1] = j;
							count_alpha2++;
						}						
					}
					else				
						rowb += b[j];
				}
			}	
//cout << "rowb= " << rowb << endl;
			// make b zero	
			if(rowb == element[1])	// == 1
			{
				if( i != 0 )
					b[leading_ones_pos[i]]=p_ONE;
				else 
				{
					if( count_1 > 0 )
					{
						mask[available_1_pos[0]]=1;
						b[available_1_pos[0]]=p_ONE;
						b[leading_ones_pos[i]]=p_ONE;
					}
					else if( count_alpha > 0 && count_alpha2 > 0 )
					{
						mask[available_alpha_pos[0]]=1;
						b[available_alpha_pos[0]]=p_ONE;
						mask[available_alpha2_pos[0]]=1;
						b[available_alpha2_pos[0]]=p_ONE;
						b[leading_ones_pos[i]]=p_ONE;
					}
				}
					
			}
			else if(rowb == element[2]) // == alpha
			{
				if(count_alpha > 0)
					mask[available_alpha_pos[0]]=1;
				else if( count_alpha2 > 0 )
				{
					mask[available_alpha2_pos[0]]=1;
					b[available_alpha2_pos[0]]=p_ONE;
					b[leading_ones_pos[i]]=p_ONE;
				}
			}
			else if ( rowb == element[3] ) // alpha^2
			{
				if(count_alpha2 > 0)
				{
					mask[available_alpha2_pos[0]]=1;
					b[available_alpha2_pos[0]]=p_ONE;			
				}
				else if( count_alpha > 0 )
				{
					mask[available_alpha_pos[0]]=1;
					b[available_alpha_pos[0]]=p_ONE;
					b[leading_ones_pos[i]]=p_ONE;
				}
			}	

			// fill the rest of positions with zero
			for(int j=leading_ones_pos[i]+1;j<n;j++)
				if( !IsZero(A[i][j]) &&  mask[j]==-1 )
					mask[j] = 1;

//cout << "mask: " << endl;
//for(int i=0;i<n;i++)
	//cout << mask[i] << " ";
//cout << endl;
//cout << "b: " << endl << b << endl;
		}
	}
	for(int i=0;i<n;i++)	// set independent variables to 1
		if(mask[i]==-1)
			b[i] = p_ONE;

	return b;
}


/**
 * given a row-reduced matrix, reduces it to reduce rwo echelon form, A_{mxn}
 */
int* reduced_row_echelon(mat_ZZ_pE& A,int p,int order)
{
	int m = A.NumRows();
	int n = A.NumCols();
	int* leading_ones_index=new int[m];
	for(int i=0;i<m;i++)
	{
		if( !IsZero(A[i]) )	// traverse non-zero rows
		{
			long j=0;
			while( j<=n && IsZero(A[i][j]) ) j++;	// find location of leading coefficient
			//  make leading coefficient a 1
			if( !IsOne(A[i][j]) )
			{
//				ZZ_pX temp_poly;
//				for(int k=0;k<A[i].length();k++)
//					SetCoeff(temp_poly,k,A[i][k]);
//				MakeMonic(temp_poly); 
				A[i]=power(A[i][j],p-2)*A[i]; // x^{p-1}=1
			}
			leading_ones_index[i]=j;	// keep track of location of leading ones in each row for sorting 
			// reduce
			for(int k=0;k<m;k++)
			{
				if(k!=i && !IsZero(A[k][j]) )
					A[k]-= A[k][j]*A[i];		// assuming leading coefficient is 1 !!!!!!!!!!!!!!!!
			}
		}
	}
	if(order)
	{
		for(int i=1;i<m;i++)
			for(int j=0;j<i;j++)
				if(leading_ones_index[i]<leading_ones_index[j])
				{
					int temp_index=leading_ones_index[i];
					vec_ZZ_pE temp_row = A[i];
					leading_ones_index[i]=leading_ones_index[j];
					A[i] = A[j];
					leading_ones_index[j]=temp_index;
					A[j] = temp_row;
				}
	}
	return leading_ones_index;
}

/**
 * Given a skew polynomial, this function recursively factors it into irreducible polynomials over the skew polynomial ring
 *
 * This function implements the Complete-Factorization algorithm in:
 * @ref "Factoring in Skew-Polynomial Rings over Finite Fields", Mark Giesbrecht
 *
 * @f skew polynomial to be factored
 * @factors array of computed factors
 *
 * @return number of factors found
 *
 **/
int SQ_CompleteFactorization(const ZZ_pEX& f,vec_ZZ_pEX& factors)
{  
  ZZ_pEX u,v;
  int iterations = 0;
  int MAX_NUM_OF_ITERATIONS=9;

  int status = SQ_Giesbretch_fac(f,to_int(ZZ_p::modulus()),(int) ZZ_pE::degree(),u,v);
  iterations++;
  while( (iterations <= MAX_NUM_OF_ITERATIONS) && (status == -1) ) // probability of failure is 8/9
    {
       status = SQ_Giesbretch_fac(f,to_int(ZZ_p::modulus()),(int) ZZ_pE::degree(),u,v);
       iterations++;
    }

  if( status == -1 )
    {
      cout << "!!! Complete Factorization Failed !!! " << endl;
      abort();
    }
  else if( status == 0 ) // f is irreducible
    {
      cout << "!!! f is irreducible !!!" << endl;
      append(factors,f);
      return 1;
    }
  else // 2 factors found, recursively find any sub-factors
    {
      ZZ_pEX g,h,R;
      h = SQ_gcrd(f,u); // h=gcrd(f,u)
      SQ_divr(f,h,g,R); // f=gh
      // factor g
      SQ_CompleteFactorization(g,factors);
      // factor h
      SQ_CompleteFactorization(h,factors);
    }
  return factors.length();
  
}

/**
 * This function is called by SQ_CompleteFactorization() to find zero divisors in a K-algebra
 *
 * @return -1 if process failed, 0 if K-algebra is a field, or else 2 on success
 * 
 **/
int SQ_Giesbretch_fac(const ZZ_pEX& f,int p,int mu,ZZ_pEX& b1,ZZ_pEX& b2,int LOG_PRINT)
{
	int n = (int) deg(f);
	int dim_of_W = n*mu;

	vec_ZZ_pE element;
	long order_of_GFq = MYgenerate_GF4(element);
	ZZ_pE theta = element[2];

	// create extension field F with primitive element \theta
	ZZ_pEX W[dim_of_W];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<mu;j++)
		{
			W[i*mu+j] = ZZ_pEX(i,power(theta,j)); // \theta^{j}x^{i}
//cout << W[i*mu+j] << endl;
		}
	}

	// Form a basis for K-algebra (K is the fixed field of F under automorphism)
	// find tranformation matrix T such that T(u)= 0 mod f for every basis element of F[x,theta]/(f)
	ZZ_pEX T_u,Q;
	mat_ZZ_p T;
	T.SetDims(dim_of_W,dim_of_W);
	for(int i=0;i<dim_of_W;i++)
	{
		SQ_divr(SQ_mulr(f,W[i]),f,Q,T_u);
		for(int j=0;j<n;j++)
		{
			ZZ_pE pE_coeff = coeff(T_u,j);	// inefficient  way of getting coefficients in NTL library
			ZZ_pX p_coeff =  rep(pE_coeff);
			for(int k=0;k<mu;k++)	// set rows of T with transformed basis vectors
				T[i][j*mu+k]= coeff(p_coeff,k);
		}
	}
//cout << "T: " << endl << T << endl;
	// find kernel of T
	mat_ZZ_p T_kernel;
	kernel(T_kernel,T);	
//cout << "kernel: " << endl << T_kernel << endl;


	// Find a non-zero zero-divisor in K-algebra
	int dim_of_Kalgebra = T_kernel.NumRows();

	// generate random a_{1} and a_{2} in K-algebra
/************************************************
*
*	Random Coefficient Generation
*************************************************/
//cout << "Enter rand vector of size " << dim_of_Kalgebra << ": ";
vec_ZZ_p rand_poly_coeff;
rand_poly_coeff.SetLength(dim_of_Kalgebra);
vec_ZZ_p x;
ZZ_pX random_poly;
random(random_poly,dim_of_Kalgebra);
for(int i=0;i<dim_of_Kalgebra;i++)			// more efficient way to assign poly coefficients to a vector ???!??!!!!
	rand_poly_coeff[i] = coeff(random_poly,i);	// 	?!?!? make sure constant coeff is non-zero ?!?!???????????!!!!!!!!!
//cin >> rand_poly_coeff;	
	mul(x,rand_poly_coeff,T_kernel);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "x: " << x << endl;
	  }
	ZZ_pEX a[3];
	for(int i=0;i < T_kernel.NumCols(); i++)
		a[0]+= x[i]*W[i];
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "a1: " << a[0] << endl;
	  }

random(random_poly,dim_of_Kalgebra);
for(int i=0;i<dim_of_Kalgebra;i++)			// more efficient way to assign poly coefficients to a vector ???!??!!!!
	rand_poly_coeff[i] = coeff(random_poly,i);	// 	?!?!? make sure constant coeff is non-zero ?!?!???????????!!!!!!!!!
//cin >> rand_poly_coeff;	
	mul(x,rand_poly_coeff,T_kernel);

	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "x: " << x << endl;
	  }

	for(int i=0;i < T_kernel.NumCols(); i++)
		a[1]+= x[i]*W[i];
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "a2: " << a[1] << endl;
	  }


SQ_divr(SQ_mulr(a[1],a[0])-SQ_mulr(a[0],a[1]),f,Q,a[2]);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "a2a1-a1a2: " << a[2] << endl;
	  }
/*************************************************/


for(int i=0;i<3;i++)
{
   ZZ_pEX b;
	// find first nu powers of b mod f where  b \in {a_{1},a_{2},a_{1}a_{2}-a_{2}a_{1}}/{0}, these powers of linearly dependent for coefficients \alpha_{0},\alpha_{1},...,\alpha_{nu}
   if(!IsZero(a[i]))
   {
	b = a[i];
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "*********************************************" << endl;
cout << "b= " << b << endl;
cout << "*********************************************" << endl;
	  }


	ZZ_pEX pows_of_b[dim_of_Kalgebra+1];
	pows_of_b[0] = SQ_power(b,0);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "b^{" << 0 << "}: " << pows_of_b[0] << endl;
	  }
	// find b^{1},b^{1},...,b^{nu}
	for(int i=1;i<=dim_of_Kalgebra;i++)
	{
		SQ_divr(SQ_mulr(b,pows_of_b[i-1]),f,Q,pows_of_b[i]); 
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "b^{" << i << "}: " << pows_of_b[i] << endl;	
	  }
	}

#if 0
      	// solve lin dep equation for \alpha_{0},\alpha_{1},...,\alpha_{nu}
	vec_ZZ_p lin_combo;
	int found_min_poly=0;
	for(int i=1;i<=dim_of_Kalgebra;i++)
	{
		mat_ZZ_pE A;
		A.SetDims(i+1,n);
		vec_vec_ZZ_pE vec_equ_system;
		for(int j=0;j<=i;j++)		// optimize: no need to make the matrix every time !!!
		{
			vec_ZZ_pE vec_equ_system_temp= VectorCopy(pows_of_b[j],n);	// use first i powers of b
			append(vec_equ_system,vec_equ_system_temp);
		}
		MakeMatrix(A,vec_equ_system);  // set up the system of equation: Ax=0
		A = transpose(A);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "A: " << endl << A << endl;
	  }
		long rank = gauss(A);
		// make each row as a polynomial monic
		

		// find solutions
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "gauss A: " << endl << A << endl;
	  }
		int* leading_ones = reduced_row_echelon(A,4,1);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "reduced A: " << endl << A << endl;
	  }
		lin_combo = solve_equ_GF2(A,leading_ones);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "solution: " << endl << lin_combo << endl;
	  }
		if( !IsZero(lin_combo) )
		{
			found_min_poly = 1;
			break;
		}
	}
#else	// temprory solution
	// try all binary combinations
	vec_ZZ_p lin_combo;
	int found_min_poly=0;
	for(int i=1;i<=dim_of_Kalgebra;i++)
	{
		mat_ZZ_pE A;
		A.SetDims(i+1,n);
		vec_vec_ZZ_pE vec_equ_system;
		for(int j=0;j<=i;j++)		// optimize: no need to make the matrix every time !!!
		{
			vec_ZZ_pE vec_equ_system_temp= VectorCopy(pows_of_b[j],n);	// use first i powers of b
			append(vec_equ_system,vec_equ_system_temp);
		}
		MakeMatrix(A,vec_equ_system);  // set up the system of equation: Ax=0
//cout << "A: " << endl << A << endl;

		
		lin_combo.SetLength(i+1);
		clear(lin_combo);
		ZZ num_of_iterations;
		power2(num_of_iterations,i);
		ZZ_p p_ONE;	
		set(p_ONE); // p_ONE =1 
		lin_combo(i+1)=p_ONE;
		ZZ counter;
		for(;counter < num_of_iterations; counter++)
		{
//cout << lin_combo << endl;
			vec_ZZ_pE product;
			product.SetLength(A.NumCols());
			for(int k=0;k<lin_combo.length();k++)
				if(!IsZero(lin_combo[k]))
					product += A[k];
			if( IsZero(product) )
			{
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << lin_combo << endl;
	  }
				found_min_poly = 1;
				break;
			}
			lin_combo = binary_inc(lin_combo);
		}	
		if( found_min_poly )
			break;	
	}
	
#endif
	// factor minimal poly
	if(found_min_poly)
	{
cout << "Found minimal polynomial ..." << endl;
		ZZ_pX min_poly;
		for(int i=0;i<lin_combo.length();i++)
			SetCoeff(min_poly,i,lin_combo[i]);	// generate corresponding polynomial from coefficients
		vec_pair_ZZ_pX_long factors;
		CanZass(factors,min_poly);
if(min_poly == factors[0].a) // if f is irreducible
{
	if( deg(min_poly) == dim_of_Kalgebra )
	{
		cout << "!!! K-algebra is a field !!!" << endl;
		return 0;
	}
}
else	// if f is reducible
{
	int num_of_DD_factors = factors.length();


	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "----------------------- min-poly factors: -------------------------" << endl;
	  }
for(int i=0;i<num_of_DD_factors;i++)
{
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
	cout << "factor " << i << ": " << factors[i] << endl; 
	  }
	ZZ_pEX fac = SQ_eval(f,factors[i].a,b);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
	cout << "f factor: " << fac << endl;
	  }
	MakeMonic(fac);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
	cout << "monic f factor: " << fac << endl;
	  }
	
ZZ_pEX Q,R;
	SQ_divr(f,fac,Q,R);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << f << " = (" << Q << ")(" << fac << ") + (" << R << ")" << endl;
	  }
}
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "----------------------------------------------------------" << endl;
	  }

	// f = gh where b1=g, b2=h
	ZZ_p p_ONE;
	set(p_ONE); // p_ONE = 1 
	ZZ_pX g=ZZ_pX(0,p_ONE);
	ZZ_pX h=ZZ_pX(0,p_ONE);

	if(num_of_DD_factors == 1 )
	{
		for(int i=0;i<factors[0].b/2;i++)
			g*= factors[0].a;		// !!???!! not taking into acount the exponents
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "g= " << g << endl;
	  }
		for(int i=factors[0].b/2;i<factors[0].b;i++)
			h*= factors[0].a;
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "h= " << h << endl;
	  }

	}
	else
	{
		for(int i=0;i<num_of_DD_factors/2;i++)
			g*= power(factors[i].a,factors[i].b);		// !!???!! not taking into acount the exponents
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "g= " << g << endl;
	  }
		for(int i=num_of_DD_factors/2;i<num_of_DD_factors;i++)
			h*= power(factors[i].a,factors[i].b);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "h= " << h << endl;
	  }
	}

	//check
	if( min_poly != (g*h) )	
		cout << "!!! ERROR: min_poly != g*h" << endl;
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "b= " << b << endl;
	  }

	b1 = SQ_eval(f,g,b);
	b2 = SQ_eval(f,h,b);
	/**********************************************
	 *                   Log
         *
         *********************************************/
	if(LOG_PRINT)
	  {
cout << "b1: " << b1 << endl;
cout << "b2: " << b2 << endl;
	  }
	return num_of_DD_factors;
}
	}
	else
		cout << "No minimal polynomial was found!!!" << endl;

   }
}
	cout << "!!! Giesbretch Factorization Failed !!!" << endl;
	return -1;
}


/**
 * Finds a right factor with of a given degree in a given skew polynomial
 *
 *
 * @ref "Factoring in Skew-Polynomial Rings over Finite Fields", Mark Giesbrecht
 *
 **/
int SQ_Giesbretch_BiFac(const ZZ_pEX& f,long s,ZZ_pEX& g,ZZ_pEX& h)
{

  // compute minimal central left multiple f_min of f
  ZZ_pX f_min1;
  ZZ_pEX f_min2;
  if ( SQ_min_central_left_multiple(f,f_min1,f_min2)!=1 )
    {
      cout << "ERROR: cannot find minimal central left multiple of " << f << endl;
      return -1;
    }
cout << "f_min= " << f_min1 << endl;
  
  // find Distinct Degree Factorization of f_min = fmin_1 ... fmin_n
  // perform Square-free decomposition
  vec_pair_ZZ_pX_long SquareFree_facs =  SquareFreeDecomp(f_min1);
  // find DDF of each Square-free factor, taking into acount the multiplicities
  vec_pair_ZZ_pX_long fmin_i;
  fmin_i.SetLength(0);

  for(int i=0;i<SquareFree_facs.length();i++)
    {
      ZZ_pX x; //x = X
      SetX(x);
      ZZ_pX h = power(x,to_int(ZZ_p::modulus()));	// h= x^p 
      vec_pair_ZZ_pX_long DDFactors =  NewDDF(SquareFree_facs[i].a,h);

      for(int j=0;j<DDFactors.length();j++)
	{
	  vec_pair_ZZ_pX_long temp;
	  temp.SetLength(1);
	  temp[0].a = power(DDFactors[j].a,to_int(SquareFree_facs[i].b));
	  temp[0].b = DDFactors[j].b;
	  append(fmin_i,temp);
	}
    }

  // find h_i = GCRD(fmin_i,f)
  vec_ZZ_pEX h_i;
  h_i.SetLength(fmin_i.length());
  for(int i=0;i<fmin_i.length();i++)
    {
      ZZ_pEX temp;     
      for(int k=0;k<=deg(fmin_i[i].a);k++) // inefficient way of converting ZZ_pX to ZZ_pEX !!!
	if( !IsZero(coeff(fmin_i[i].a,k)) )
	  SetCoeff(temp,k);
	  
      h_i[i]= SQ_gcrd(temp,f);
    }

  // Factor each h_i completely in skew polynomial ring as h_i = h_i_1 ... h_i_ei
  int num_of_h_polys = h_i.length();
  vec_ZZ_pEX h_i_facs[num_of_h_polys];
  int e_i[num_of_h_polys],i_list[num_of_h_polys];

  for(int i=0;i<num_of_h_polys;i++)
    { 
      SQ_CompleteFactorization(h_i[i],h_i_facs[i]);
      e_i[i]= h_i_facs[i].length();
      i_list[i] = (int) deg(h_i_facs[i][0]);
    }

//cout << "h len= " << num_of_h_polys << endl;


  // determine the set d_1,...,d_n
  int* d_i; 
  d_i = new int [num_of_h_polys];
  //reset
  for( int ppp=0;ppp<num_of_h_polys;ppp++)
	d_i[ppp]=0;
  

  int status =  int_decomposition(i_list,e_i,s,num_of_h_polys,d_i);

/*
int sum=0;
for(int jjj=0;jjj<num_of_h_polys;jjj++)
{
 cout  <<  i_list[jjj] << "*" << d_i[jjj] << " + ";
 sum += i_list[jjj]*d_i[jjj];
}
cout << endl << "sum = " << sum << endl;
cout << "s= " << s << endl;
if(sum != s )
cout << "!!! ERROR: sum != s !!!" << endl;
*/

  if(status)
    {
      // find sub_h_i's
      vec_ZZ_pEX sub_h_i;
      sub_h_i.SetLength(num_of_h_polys);
      for(int mm=0;mm<num_of_h_polys;mm++)
	{
	  SetCoeff(sub_h_i[mm],0); // sub_h_i = 1

	  for(int nn=0;nn<d_i[mm];nn++)
	    {
//cout << "h_i_facs= " << h_i_facs[mm][e_i[mm]-nn-1] << endl;
	      sub_h_i[mm] = SQ_mulr(h_i_facs[mm][e_i[mm]-nn-1],sub_h_i[mm]);
	    }
	}

      // h = lclm(sub_h_1, ... , sub_h_n)
      if( num_of_h_polys == 1 )
	h =sub_h_i[0];
      else
	{
	  h = SQ_lclm(sub_h_i[0],sub_h_i[1]);
	  for(int mm=2;mm<num_of_h_polys;mm++)
	    h = SQ_lclm(h,sub_h_i[mm]);
	}
cout << "deg(h)= " << deg(h) << endl;

      // find g such that f=gh
      ZZ_pEX RRem;
      SQ_divr(f,h,g,RRem);

cout << "h= " << h << endl;
cout << "g= " <<  g << endl;

      if( !IsZero(RRem) )
	{
	  cout << "!!! ERROR: RRem != 0: Bi-factorization of f is incorrect !!!" << endl;
	  return -1;
	}
      if( deg(h) != s )
	{
	  cout << "!!! deg(h)!= s !!! " << endl;
	  return -1;
	}

      return 0;
    }

  cout << "!!! f has no right factor of degree " << s << " !!! " << endl;
  return -1;


}

/**
 * Finds minimal central left multiple of a given skew polynomial
 *
 * Minimal Central Left Multiple f_min in K[x] (fixed field of F) of a polynomial f in skew ring F[x,\sigma]  is the polynomial of the minimum degree such that f_min = 0 (mod f) 
 *
 **/
int SQ_min_central_left_multiple(const ZZ_pEX& f,ZZ_pX& f_min1,ZZ_pEX& f_min2)
{
  long n = (long) deg(f);
  int mu = (int) ZZ_pE::degree();

  vec_ZZ_pEX R;
  R.SetLength(n*mu+1);
  set(R[0]); // R[0] = 1

  ZZ_pEX x_mu;
  SetCoeff(x_mu,mu);

  mat_ZZ_p B;
  B.SetDims(n*mu+1,n*mu);
  for(int j=0;j<n;j++)
    {
      ZZ_pE pE_coeff = coeff(R[0],j);	// inefficient  way of getting coefficients in NTL library
      ZZ_pX p_coeff =  rep(pE_coeff);
      for(int k=0;k<mu;k++)	// set rows of T with transformed basis vectors
	B[0][j*mu+k]= coeff(p_coeff,k);
    }

  for(long i=1;i<=n*mu;i++)
    {
      ZZ_pEX Q;
      SQ_divr(SQ_mulr(x_mu,R[i-1],&mu,1),f,Q,R[i]);
 // cout << "R[" << i << "]= " << R[i] << endl;
      for(int j=0;j<n;j++)
	{
	  ZZ_pE pE_coeff = coeff(R[i],j);	// inefficient  way of getting coefficients in NTL library
	  ZZ_pX p_coeff =  rep(pE_coeff);
	  for(int k=0;k<mu;k++)	// set rows of T with transformed basis vectors
	    B[i][j*mu+k]= coeff(p_coeff,k);
	}

    }


  //cout << "B= " << endl << B << endl;

// temprory solution
// try all binary combinations
	vec_ZZ_p lin_combo;
	int found_min_poly=0;
		
	lin_combo.SetLength(n*mu);
	clear(lin_combo);

	ZZ num_of_iterations;
	power2(num_of_iterations,n*mu);

	ZZ_p p_ONE;	
	set(p_ONE); // p_ONE =1 
	lin_combo(1)=p_ONE;

	ZZ counter;
	for(;counter < num_of_iterations-2; counter++)
	{
//cout << lin_combo << endl;
		vec_ZZ_p product;
		product.SetLength(B.NumCols());
		for(int k=0;k<lin_combo.length();k++)
			if(!IsZero(lin_combo[k]))
				product += B[k];
		if( IsZero(product) )
		{
//cout << "lin_combo= " << lin_combo << endl;
			found_min_poly = 1;
			break;
		}
		lin_combo = binary_inc(lin_combo);	
	}

	if(found_min_poly)
	  {
	    for(int k=0;k<lin_combo.length();k++)
	      if(!IsZero(lin_combo[k]))
		{
		  SetCoeff(f_min1,k*mu);
		  SetCoeff(f_min2,k*mu);
		}
	  }

  return found_min_poly;
}
